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ABSTRACT 



We carry out three-dimensional magnetohydrodynamical simulations of the 
magnetorotational (Balbus-Hawley) instability in weakly-ionized plasmas. We 
adopt a formulation in which the ions and neutrals each are treated as separate 
fluids coupled only through a collisional drag term. Ionization and recombination 
processes are not considered. The linear stability of the ion-neutral system 
has been previously considered by Blaes & Balbus (1994). Here we extend 
their results into the nonlinear regime by computing the evolution of Keplerian 
angular momentum distribution in the local shearing box approximation. We 
find significant turbulence and angular momentum transport when the collisional 
frequency is on order 100 times the orbital frequency Q. At higher collision rates, 
the two-fluid system studied here behaves much like the fully ionized systems 
studied previously. At lower collision rates the evolution of the instability is 
determined primarily by the properties of the ions, with the neutrals acting 
as a sink for the turbulence. Since in this regime saturation occurs when the 
magnetic field is superthermal with respect to the ion pressure, we find the 
amplitude of the magnetic energy and the corresponding angular momentum 
transport rate is proportional to the ion density. Our calculations show the 
ions and neutrals are essentially decoupled when the collision frequency is less 
than 0.01^2; in this case the ion fluid behaves as in the single fluid simulations 
and the neutrals remain quiescent. We find that purely toroidal initial magnetic 
field configurations are unstable to the magnetorotational instability across the 
range of coupling frequencies. 



Subject headings: accretion, accretion disks-protostellar disks-instabilities-MHD 
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1. Introduction 

The key to understanding accretion disk dynamics lies with the angular momentum 
transport mechanism. Since the molecular viscosity of disks is very low, some form of 
"anomalous viscosity" must be present. Although the precise nature of this anomalous 
viscosity has long been elusive, the discovery that differentially rotating systems are 
magnetohydrodynamically (MHD) unstable (Balbus & Hawley 1991) has led to the 
conclusion that fully ionized disks must be MHD turbulent. Because this magnetorotational 
instability is caused by angular momentum transport, the resulting turbulence has precisely 
the right character to transport angular momentum outwards (Balbus, Hawley & Stone 
1996, hereafter BHS96), as required for disks to accrete. Since the MHD instability plays 
such a fundamental role in disks, the question naturally arises, what is its behavior when 
the plasma is not fully ionized? Protostellar and protoplanetary disks as well as other 
molecular disks are the venues where this question is particularly significant. 

Although it is tempting to assume that in the absence of MHD turbulence, purely 
hydrodynamical turbulence will rise to the task of transporting angular momentum, this 
now seems highly unlikely. BHS96 demonstrated, through a combination of analysis and 
simulation, that differentially rotating systems are both linearly and nonlinearly locally 
stable to hydrodynamic perturbations so long as the standard Rayleigh criterion is satisfied. 
Even if initiated "by hand," hydrodynamic turbulence is not self-sustaining. Outward 
transport of angular momentum through a net Reynolds stress requires a specific average 
correlation between the radial velocity and the angular momentum fluctuations in the 
turbulent flow. Of fundamental importance is the interaction of these velocity fluctuations 
with the background mean flow. In differentially rotating systems the source of free energy 
is the angular velocity gradient. Angular momentum fluctuations, on the other hand, act 
to reduce the background angular momentum gradient, which has the opposite sign from 
the angular velocity gradient. A positive value of the Reynolds stress, required to tap 
into the free energy of the system, acts as a sink term for the evolution of the angular 
velocity fluctuations that make up the Reynolds stress itself. Thus the turbulence is not 
self-sustaining. 

The results of BHS96 further imply that enhanced angular momentum transport is not 
the necessary outcome of turbulence. Because transport requires a high degree of correlation 
between the radial and azimuthal velocity fluctuations, turbulence that does not have its 
origin in the mean differentially rotating flow is unlikely to be an efficient source of angular 
momentum transport. This point is emphasized by Stone & Balbus (1996), whose numerical 
simulation showed that while vertical convection can generate turbulence, the resulting 
net radial angular momentum transport was inward at a low rate (so long as the equator 
was kept hot by the boundary conditions of the simulation). Similar results were found by 
Cabot (1996). Simulations by Ryu & Goodman (1994) of the action of a parametric tidal 
instability on a disk show the generation of turbulence, but without internal transport. 
These studies provide compelling evidence that the idea of purely hydrodynamic turbulence 
as an angular momentum transport mechanism in protostellar disks should be discarded. 
If turbulence is to be driven by the differential rotation, it must be MHD turbulence. If 
hydrodynamic turbulence is present, it must be driven by some source other than the 
background shear and generally will not transport angular momentum. Such considerations 
are particularly important if turbulence is important in protoplanetary disks for chondrule 
and planetesimal formation (e.g., Cuzzi, Dobrovolskis & Hogan 1996). 

It now appears that the number of different mechanisms for transporting angular 
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momentum in disks is very limited. Magnetic fields and spiral wave disturbances 
(nonaxisymmetric distortions) are the only viable general mechanisms. In the latter 
category, and within the context of protostellar disks, nonaxisymmetric gravitational 
instabilities may be important. Global wave mechanisms, however, are not normally 
associated with the generation of local turbulence. The dynamics of a disk controlled 
by such global waves will be quite different from a model based upon the usual local 
a prescription for viscosity, which is itself based upon the ansatz of transport by local 
turbulent stresses. 

For protostellar disks we are forced to reckon with questions of the effectiveness of 
magnetic field coupling. There is, of course, nothing about this issue that is unique to the 
magnetorotational instability. If the magnetic field is not well coupled to the fiuid, then all 
potential MHD processes in the disk will be affected. It is difficult to imagine, for example, 
a scenario in which a weak field in a disk is stablized by low ion-neutral coupling, yet 
remains involved in, say, a kinematic dynamo. If a protostellar disk or rotating molecular 
cloud is sufficiently well-coupled to a weak magnetic field for the field to be important in 
any way, then the magnetorotational instability is an important dynamical factor. 

Under what circumstances, then, will the instability operate and produce MHD 
turbulence in protoplanetary disks where the ionization fraction is quite low? The first 
step towards answering this question was taken by Blaes & Balbus (1994, hereafter BB), 
who performed an axisymmetric linear stability analysis for a vertical (and toroidal) field 
in a weakly ionized plasma for a number of limits. Their major finding is that the linear 
magnetorotational instability is present so long as the ion-neutral collision frequency 
exceeds the local epicyclic frequency. This condition will be satisfied even for very small 
ionization fractions. They also found that azimuthal fields can reduce the observed growth 
rates, although such fields do not eliminate the instability. 

While a linear analysis can indicate the presence or absence of the instability, its 
effectiveness as a transport mechanism must be determined by its nonlinear evolution. The 
first numerical study of the two-fluid magnetorotational instability was carried out by Mac 
Low et al. (1995). They assume ionization-rccombination equilibrium and consider the 
low-ionization limit, neglecting the ion pressure and inertia. This reduces the problem to 
a single-fluid (neutral) plus a diffusion term in the induction equation (ambipolar diffusion 
limit). Using the ZEUS code they carried out a series of two-dimensional simulations of the 
vertical flux tube problem (as in Hawley & Balbus 1991) for various ion-neutral coupling 
strengths. Although these simulations did not follow the evolution much beyond the linear 
stage, the results were in agreement with the stability analysis of BB in the appropriate 
limit. Essentially, so long as there are unstable wavelengths available, and the coupling 
between ions and neutrals is sufficiently strong, the instability behaves much like the single 
fluid case. The instability ceases to operate when the ambipolar diffusion rate becomes 
comparable to the growth rate of the instability, i.e., the fleld diffusion time is < where 
is the disk orbital frequency. 

In another study, Brandenburg et al. (1995) investigated the ambipolar diffusion limit 
in three-dimensional simulations of a local, vertically stratified disk. They considered a case 
where the ambipolar diffusion time was long compared to the orbital time. They found that 
in this limit (i.e., ambipolar diffusion sufficiently small) the instability remains effective, 
and continues to generate self-sustained turbulence that transports angular momentum 
outward, albeit at a rate slightly reduced from the fully-coupled case. In another simulation, 
the diffusion time was set comparable to and the turbulence decayed. 
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These first results, while important, are only a beginning. To date, all the numerical 
studies have considered the behavior of the partially ionized system in the limit where the 
inertia of the ions can be completely neglected, and where the ion density is everywhere a 
fixed power-law function of the neutral density. In this paper we will approach the problem 
using a genuine two-fluid, ion-neutral evolution to examine effects where the ions are free 
to move relative to the neutrals. We will investigate the transition from the well-coupled 
regime, through critical coupling where the collision frequency is comparable to the epicyclic 
frequency, and down to the fully uncoupled limit. This parameter study should more 
clearly define the physical conditions for which full MHD turbulence and accompanying 
angular momentum transport can be expected, and those for which the bulk of the system 
dynamics are essentially hydrodynamic. The full range of conditions is likely to be of 
importance somewhere within protostellar or protoplanetary disks. In some regions of these 
disk systems, the ionization fractions may be quite small, but in other regions, such as 
near the forming protostar, the temperatures will be high, and the gas will be nearly fully 
ionized. In between, there will be a transition region. Determining the size of this transition 
region, and delineating its properties, will depend on obtaining a better understanding of 
the nonlinear behavior of the ion-neutral system in various regimes. 

The plan of the paper is as follows. In §2 we will consider the equations and 
the numerical techniques used to solve them. Although the collision term is handled 
semi-implicitly, the use of explicit finite-differencing for the remainder of the system makes 
the code Courant-limited. As a test problem, we compare numerical growth rates with 
analytic values from BB. In §3 we present the results of an extensive ensemble of simulations, 
covering a range of ionization fractions and coupling frequencies. Because of the Courant 
limit, this study is limited to relatively large ionization fractions / = Pi/{pi + p„)- In the 
linear limit, however, the growth rates for small ionization fractions are relatively unaffected 
by decreasing /, if the ratio of the coupling frequency to orbital frequency is held constant. 
Thus, we expect that physical insights gained by this study should extend even into the 
small / regime. The implications of the simulations will be summarized and discussed in §4. 



2. Two-Fluid Algorithm 



2.1. Equations and Numerics 



In these simulations we will be studying the simple two-fluid system as described 
in §3 of BB, consisting of separate compressible ion and neutral fluids, coupled only by 
collisions. The ion fluid is assumed to be perfectly conducting. We include no ionization 
or recombination; ions and neutrals are separately conserved. Wc have examined both 
isothermal systems, and ones where the two fluids can have different temperatures and no 
thermal coupling is included. 

As in previous work (Hawley, Gammie & Balbus 1995, hereafter HGB95) we restrict our 
attention to the local Hill (1878) system representation of a disk. This model incorporates 
Coriolis and tidal forces, and is constructed by expanding the equations of motion in 
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cylindrical coordinates (R, 0, z) to first order around a fiducial radius Ro using a set of local 
Cartesian coordinates x = (x, y, z) = {R — Ro, Ro[(f> — ^t], z). The rotation law is assumed 
to have the form Vt ~ for a Keplerian disk q = 3/2. The resulting equations for the 
ions, neutrals, and magnetic field are 

+ V • (av,) = 0, (1) 



dt 

dpn 

dt 



+ V • (pnV„) = 0, (2) 



--^ + • Vvi = V {Pi + ^) + ^ " ^ - 20 X + 2g02xx + 7p„(v„ - v^), (3) 

^ + v„ • Vv„ = VP„ - 2J2 X v„ + 2g02^x - 7Pi(v„ - v^), (4) 

ut Pfi 

— = V X (v, X B), (5) 

where the terms have their usual meaning, the subscripts n and i refer to the neutrals and 
ions respectively, and 7 is the drag coefficient due to collisions between the two fluids. In 
addition, there is an equation of state (we use either an adiabatic or an isothermal equation 
of state) for the neutrals and for the ions. In these equations we neglect vertical gravity and 
assume that the ion component is a perfectly conducting fluid. While there are physically 
important regimes where resistivity will be important (see, e.g., the discussion in BB), we 
make this assumption to simplify the parameters of the present study. The results of this 
paper will correspond to the most favorable conditions for the MHD instability since finite 
resistivity would be expected to reduce the levels of the resulting turbulence. 

Equations (l)-(5) amount to separate hydrodynamical (neutrals) and MHD (ions) 

systems, coupled by the collisional drag term. We use the hydrodynamical algorithms 
described by Stone & Norman (1992a) to evolve the equations for the neutrals (eqs. [1-2] 
and [4]), and the MOCCT algorithm (Stone & Norman 1992b; Hawley & Stone 1995) 
to evolve the ions (eqs. [3] and [5]). The collisional drag term 7 in both equations (3) 
and (4) is operator-split and updated using fully implicit backward Euler differencing 
(a complete description of the resulting difference equations is given in Stone 1997). 
This avoids the Courant stability criterion associated with diffusion of the magnetic field 
that restricts explicit differencing formalisms of the drag term. The resulting hybrid 
hydrodynamical-MHD code has been tested with steady C-type shock solutions and 
comparison to the analytic growth rate of the Wardle instability (Stone 1997). We describe 
a comparison of the numerical growth rate of the magnetorotational instability in a partially 
ionized disk with the analytic value from BB in §2.2 below. 

The computations are done in the periodic three-dimensional (3D) shearing box 
(HGB95). In this system the computational domain is a rectangular box with sides L^jLy, 
and Lz, and it is assumed to be surrounded by identical boxes that are strictly periodic at 
t = 0. A large-scale continuous hnear shear flow is present across all the boxes. At later 
times the computational box remains periodic in y and z, while the radial (x) boundary 
condition is determined by the location of the neighboring boxes as they move relative to 
one another due to the background shear. We use box dimensions of L^. = 1, Lj^ = 27r, and 
L2 = 1, as in previous single- fluid simulations (HGB95). 
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The initial equilibrium system has a Keplerian rotation law {q = 3/2); the angular 
speed at the center of the shearing box is fl. The initial state has constant densities 
and pressures (ion and neutral). For most experiments we use an isothermal equation of 

state; in a few simulations we use of an adiabatic equation of state with P oc p^^^. The 
mean molecular weight of the ions and the neutrals can be specified separately; for most 
simulations we adopt the weights used in BB, m„ = 2.33mH and = 30m//. These 
numbers come from assuming that the neutrals are hydrogen and helium molecules while 
the ions correspond to trace alkali species (Draine, Roberge, & Dalgarno 1983). Assuming 
thermal equilibrium, the ratio of the initial sound speeds in the neutrals and ions is equal 
to the square root of the ratio of the mean molecular weights. 

One of the difficulties of this study is the wide range of possible physical parameters 
that could be investigated. Among these are the equation of state, neutral to ion mass 
ratio, initial field strength and orientation, the initial ionization fraction / = Pi/{pi + Pn), 
and the coUisional frequency, Pij/fl. Our primary focus will be on the role of the coUisional 
drag term 7. In the limit of 7 — 0, the ion and neutral fluids decouple; the ions evolve 
as in the magnetized single-fluid case, while the neutrals evolve hydrodynamically. In the 
other extreme, 7 — > 00, the system reduces to a single well-coupled fluid. In this study, 
wc will examine a range of coUisional frequencies on either side of the critical frequency, 
■ypi ~ fl, by treating 7 as a free parameter and varying it for given values of ionization 
fraction and Q. The Courant condition from the ion Alfven speed limits the simulations to 
moderate values of the ionization fraction / > 0.001. As BB have shown, it is the ratio of 
the coUisional frequency to the orbital frequency that determines the linear behavior of the 
instability. Assuming this holds even into the nonlinear regime, our results should be able 
to set important limits even for significantly smaller, and hence more realistic, ionization 
fractions. 



2.2. Two-Dimensional Test Simulation 



As a simple verification test of the two-fiuid code, we compare observed growth rates 
in the linear regime of the instability with analytic values from BB. (Similar tests were 
performed for the fully-coupled, single-fiuid system by Hawley and Balbus 1992.) We 
choose an initial magnetic field strength corresponding to Pi = Pion/Pmag = 2, run several 
models with / = 0.1 and different values 7, and compare the observed growth rate in the 
perturbed magnetic radial field with the value given by the dispersion relation (BB 
eq. [25]) for the largest vertical wavelength permitted by the computational domain. The 
particular initial magnetic field strength used makes the (analytically determined) fastest 
growing wavelength close to the domain size = 1. We use two resolutions, 63 x 63 and 
31 X 31 grid zones; these are the {x, z) resolutions that will be used in the 3D simulations to 
follow. The numerical growth rates, in units of Q, arc plotted on top of an analytic growth 
rate curve obtained from equation (24) of BB (Fig. 1). As can be seen from the figure, 
the two-fiuid code reproduces the analytic linear growth rates for the range of collision 
frequencies considered. There is not much difference between the high and low resohition 
simulations; the long wavelength modes considered here are adequately resolved by either 
grid. 
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3. 



Results 



3.1. A Review of Single Fluid Results 



In this paper our primary focus will be on an initially uniform vertical magnetic field, 
Bz- The single- fluid, fully-ionized version of this problem was investigated by HGB95, and 
we begin with a brief summary of their results. 

The maximum growth rate of the single-fluid system is 0.75Q. HGB95 found this 
growth rate for wavenumbers kv^z/^ ~ 1 in all simulations where the fastest growing 
wavelength is adequately resolved on the grid. The fastest growing mode is axisymmetric 
and leads to flow along radial "channels." Outward flowing channels contain excess angular 
momentum, while inward flowing channels have less than the Keplerian value. When the 
fastest growing mode has a vertical wavelength greater than or equal to the radial dimension 
of the box, this axisymmetric mode continues to grow exponentially well into the strong 
field regime. However, for smaller vertical wavelengths (larger vertical wavenumbers), the 
channels break down into turbulence. These numerical results can be understood through 
the analysis of Goodman & Xu (1994) who found that, in the local limit and for a uniform 
vertical field, the exponentially-growing linear solution is also an exact nonlinear solution. 
Goodman & Xu further carried out a stability analysis of the streaming channel solution 
and found that it was subject to a number of instabilities. The most important of these 
parasitic instabilities is a magnetized Kelvin-Helmholtz mode that is present for radial 
wavelengths larger than the channel solution's vertical wavelength. Thus, under most 
circumstances, the channel solution is unstable, and the general outcome of the instability 
is turbulence. 

The resulting turbulence is nonisotropic; correlated fiuctuations in the magnetic and 
velocity fields transfer angular momentum outward through the action of the Maxwell 
(magnetic) and the Reynolds stresses. Accretion disk models often make use of the Shakura 
& Sunyaev (1973) parameterization for the stress Wr^, scaling it with the total pressure, 
Wr((, — OiP. In the MHD turbulence, however, the level of net transport depends primarily 
upon the magnetic field levels in the saturated turbulent state. The stress is proportional 
to the magnetic pressure with a ratio of the stress to Pmag of about 0.6. We refer to this 
ratio as ctmag- 

The complete ensemble of HGB95 vertical field simulations provides an empirical 
relation for the average magnetic field energy in the turbulence: it is proportional to the 
product of the fastest growing initial wavelength and the maximum possible unstable 
wavelength permitted in the computational domain. From HGB95, we have 



where Ac is the fastest growing wavelength for the initial magnetic field 
(Ac = 27r[16/15]i/2^A/^^)- A best fit to the HGB95 simulations gives a constant of 
proportionality of 1.2. In essence, the final energy is the geometric average of the energy of 
the net vertical field piercing the box (which, because of the periodic boundary conditions 
does not change as a result of the evolution) and the energy of the largest magnetic field 




(6) 
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that would be unstable given the dimension L^. The box dimension was set equal to 
Cg/O. to represent a disk scale height, although vertical gravity was not included in these 
simulations. The growth of the instability and its saturation level are unaffected by the 
hydrodynamic pressure in the box; the magnetic field energy remains subthermal (since the 
dimensions of the computational domain were chosen so that the sound speed was equal to 
L^O.) The saturation amplitude is also unaffected by the presence of a subthermal toroidal 
field except to the extent that such a toroidal field and its instabihties add to the total 
magnetic energy. 



3.2. Two-fiuid system: uniform fields 



We now turn to the two-fluid system and begin with a simple initial field geometry, a 
uniform Bz field; in subsequent sections we consider the evolution of B^ fields with zero net 
fiux, uniform By fields, and mixtures of these cases. 

In the two-fluid systems there are two identifiable limits in which the single-fluid 
results should apply: 7 = and 7 00. When 7 = the ions are completely decoupled 
from the neutrals and their evolution will be as a single fluid with the ion density and ion 
Alfven speed determining the final magnetic energy levels. In the other limit, the system 
will again behave as a single fluid, but with the iMal density determining the Alfvcn speed 
and saturation energies. BB point out that the linear properties of the two-fluid system 
can be understood in terms of the "effective" density of the system, which varies between 
Pi and Pi + pn depending on the strength of the coupling. Since \c oc va = B/Airp^^"^, the 
scaling relation (6) implies that the ratio of the magnetic energy in the decoupled (7 = 0) 

limit to the magnetic energy in the fully coupled (7 = 00) limit will be Z^^^, for fixed initial 
magnetic field strength. We will explore the properties of the transition region between 
these two limits by varying the drag coefficient 7, the ionization fraction /, and other 
variables such as the equation of state and initial field strength in a number of different 
models. 



3.2.1. A fiducial run 



We begin our discussion with a baseline simulation against which to compare other 
models. This fiducial run (Z17) is a two-fiuid shearing box with 31 x 63 x 31 resolution, 
the standard survey resolution. (The complete list of vertical field simulations is provided 
in Table 1.) We set the neutral density p„ = 1, and the neutral sound speed Cg = Lz^, 
with Lz = 1 and Q = 10~^. The ion fraction is set to / = 0.1. This, together with a 
the drag coefficient 7 = 0.01, yields a coupling frequency of 'jpi/fl — 1.11, just above the 
critical value. We use an isothermal equation of state, and assume a ratio of ion to neutral 
mass nii/mn of 30/2.33 — 12.9. The ion sound speed is equal to Cgi — {rrii/mnY/'^. The 
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assumption that the ions and neutrals have the same temperature estabhshes the relative 
pressures. The initial field is a uniform field with strength /3j = 80, which corresponds to 
a value of kvAi/^ = O.28L2/A. From the linear dispersion relation [eq. (24) of BB] we find 
that the unstable growth rates increase with increasing wavenumber, from a rate of 0.15J1 
for X = Lz to OASfl for A = I/2/4. Smaller wavelengths are also unstable with comparable 
growth rates, up to the limit of A = O.OGSLz, but such wavelengths are not well resolved. 

The uniform initial conditions are perturbed with small, random fluctuations in the 
ion and neutral pressures. The magnetic fleld energy grows exponentially for the first 
three orbits before leveling out after orbit 4 (Fig. 2). The numerical growth rate of the 
instability is computed from the slope of the radial field energy density curve. The radial 
field perturbation grows at a rate of 0.270, equal to the analytic value for A = Lz/2. 
At 4.6 orbits plots of variables such as the magnetic field show significant A = L2/4 
structure. At this point, the system is largely axisymmetric, having the appearance of the 
channel solutions previously seen in pure B^ single fluid simulations (Hawley & Balbus 
1992; HGB95). 

At orbit 15 the toroidal fleld is dominated by a nearly axisymmetric \ = structure. 
The ions are greatly affected by the fleld since the toroidal fleld is now superthermal with 
Pi — 0.6. Filaments of high density ions are sandwiched between these regions of strong 
toroidal fleld (Fig. 3). The maximum ion density in the fllaments is 2.5 times the average 
value. Since the neutral density has a much smaller relative variation, the ionization fraction 
in the fllaments has gone up by a comparable amount. The neutral density plot shows much 
less small scale structure; the neutrals are dominated by trailing m = 1 pressure waves 
that are nearly independent of z. These have a maximum 5Pn/Pn of about 10%. Thus 
considerable separation of the ions from the neutrals occurs when the coupling frequency is 
comparable to the orbital frequency. 

Figure 2 shows signiflcant magnetic energy fluctuations after saturation. Representative 
values are obtained by time and space averages; some of these values are listed in Table 
2. The time-averaged magnetic fleld energy after saturation is By/Sn — 0.011{LzO,y; the 

majority of the energy is in the toroidal fleld with B^/B^ = 0.048, and B^/B^ = 0.035. 
Since the toroidal fleld grows by both the action of the instability and the background shear 
flow, the amplitude of the poloidal flelds provides the best measure of the turbulence. In 
this simulation the poloidal flelds are considerably smaller than in single fluid runs. For 
example, in the flducial run of HGB95 B^/By = 0.32, and Bl/By = 0.10. Interestingly, the 

average magnetic fleld energy in the two-fluid system lies below both of the two single-fluid 
limits given by equation (6), either the one appropriate for the ions alone, or the one 
appropriate to the total density of ions plus neutrals. This implies that near the critical 
coupling frequency the neutrals primarily act as a drag to reduce the vigor of the fleld 
ampliflcation, the MHD turbulence and the resulting transport. 

The time-averaged perturbed kinetic energies in the ions are \{pv1)i = 6.8 x 10~'^{LzQ)^, 
\{p6vl)i = 3.0 X \Q-\LzQ.f, \{(yvl)i = 7.0 x IQ-^iL^nf. The evolution of the ion and 
neutral kinetic energy is depicted in Fig. 4. The volume-averaged kinetic energy of the 
neutrals evolves closely with that of the ions, except that it is larger by nearly the ratio 
pn/ pi- The ion velocities are slightly larger on average. 

As in the single-fluid case, angular momentum is transported outward by the Reynolds 
stress (both neutral and ion) and by the Maxwell stress. The Maxwell term is the largest 
stress component with —{B^By/An) — 0.0034(L2Q)^; the neutral and ion Reynolds stresses 
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are 2.4 x 10~^ and 3.5 x 10"'^{Lzfl)'^ respectively. The ion Reynolds stress is proportionally 
larger than that of the neutrals, once the neutral to ion density ratio is factored out. 
In these units the sum of the stresses is equal to the Shakura-Sunyaev a value. Single 
fluid simulations find that while the ratio of total stress to the total pressure can vary 
substantially from one simulation to another, the ratio of the stress to the magnetic pressure 
is more nearly constant. The direct dependence of the stress on the magnetic pressure 
reflects both the fundamental role that the magnetic field plays in driving both the Maxwell 
and Reynolds stresses, as well as the high degree of correlation between the radial and 
toroidal field components. Clearly this remains true for ion-neutral systems. In the fiducial 
run ttmag = 0.61, a value comparable to that seen in the single-fluid ^-field simulations. 
Despite its lower relative value compared with the single-fluid simulation, the radial fleld 

remains highly correlated with the toroidal. The ratio (BxBy) / ^{B"^.) (By) is nearly 1 

during the linear growth phase and has an average value of 0.7 after saturation. 

The existence of unstable wavenumbers with large values, and the relatively low 
growth rate observed in the magnetic field, suggests that resolution may be infiuencing the 
simulation. Resolution can be especially important for ion-neutral simulations, because 
they typically have a wider range of rapidly growing wavelengths and physically important 
lengthscalcs compared to the single fluid system. For the parameters of the flducial run, 
this is certainly true, and more grid zones should be required to resolve the linear growth 
phase adequately compared to the one-fluid system. Run LZl doubles the resolution of the 
fiducial run to 63 x 127 x 63 grid zones while maintaining all other initial parameters. The 
observed growth rate for the magnetic field in the high resolution run is 0.35f2 (see Fig. 
2), 30% larger than in the lower resolution run. The most obvious mode seen in density 
plots during the linear growth phase has L^/X = 6. Although ultimately both the low 
and high resolution simulations are dominated by the largest scales, resolution remains an 
important factor throughout the simulation. The high resolution simulation has 1.8 times 
as much magnetic energy as the fiducial run, when compared over the same time interval 
after saturation. The total stress is larger by a factor of 1.4. 



3.2.2. Effect of ionization fraction 

In the next set of simulations, we investigate the effect of the ionization factor by 
comparing the fiducial run to model Z20 with / = 0.01 and 7 = 0.1, and model Z21 with 
/ — 0.001 and 7 = 1. We retain the same initial Alfven speed vm = O.OAA{Lz^l) as used in 
the fiducial run. Because we hold the Alfven speed constant, the initial magnetic field 
decreases as f^^^. Table 2 lists time and space averaged values after saturation for this set 
of runs. 

Although the parameter kvAi/^ = O.28L2/A is the same from one run to the next, 
the linear growth rates for the resolvable vertical wavelengths decrease with / (BB). For 
X = Lz the linear growth rates are 0.15, 0.047, and 0.015^7; for A = Lz/2 the rates are 0.27, 
0.093, and 0.03f2. The measured growth rates for B^ in the simulations are 0.27, 0.072, and 
0.016Q for / = 0.1, 0.01, and 0.001 respectively. 

Despite the reduced growth rates, the total magnetic field amplification is comparable 
in all three runs, about two orders of magnitude in energy. The time-averaged magnetic 
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energy after saturation in these runs is 1.2 x 10 ^, 1.1 x 10 ^, and 2.0 x 10 ^(L^fi)^. These 
correspond to similar saturation amphtudes in the ion Alfven speed, ~ 0.1 (L^l^)^ (Fig. 
5). The majority of the magnetic energy is in the toroidal component. Since the saturated 
magnetic energy is proportional to /, the total stress is similarly proportional to / as is 
the Shakura-Sunyaev a parameter, defined with respect to the total pressure (neutral plus 
ion pressure). More significantly, the relative strength of the radial field compared to the 
toroidal declines with /. This means that the average ctmag ^-Iso declines with /; ctmag = 0.6, 
0.2, and 0.1 in these three runs. 

As in the fiducial run, the ions arc squeezed into vertically thin sheets lying between 
regions of strong toroidal field; the lower the initial ionization fraction, the greater the 
compression. At a representative late time in runs Z17, Z20, and Z21, the maximum ion 
densities are 2.5, 2.6 and 4.2 times the initial value. The neutrals look similar in all three 
runs, again dominated by nearly 2;-independent trailing pressure waves. 

Because the three models saturate at comparable ion Alfven speeds, we conclude that 
the final state depends most strongly upon the ion density, that is, the effective inertia for 
the instability remains close to pi. At coupling frequencies ~ 1 the role of the neutrals is 
mainly that of a sink for the turbulence. This is consistent with the observation that the 
magnetic energies in these simulations are smaller than predicted from equation (6) for a 
single-fluid with density equal to pi. 

While these results imply an increasingly small angular momentum transport for 
increasingly small ionization factors, this is partly because we have chosen to initialize 
this ensemble of simulations with constant Alfven speed. An alternative is to begin with 
constant magnetic energy density. For example, a model with / = 0.001 is still unstable 
even with the field is as strong as in the fiducial / = 0.1 simulation. With this strength of 
magnetic field and with / — 0.001, we have Pi — 0.67 initially, and kvAi/^ — S.SSL^/A. The 
analytic growth rate is 0.35 for the X — Lz mode. When such a run (Z23) is performed, 
however, the observed growth rate is only 0.150 compared with 0.270 in the fiducial 
run. The field is sufficiently strong (relative to the ion pressure), that compressibifity is 
important, and this reduces the growth rate below the weak-field, incompressible value 
obtained from equation (24) of BB. 

The saturated magnetic field energy in Z23 is lower than in the fiducial run, with 
B'^/Stt = 1.3 X 10^^{Lz^)'^. This is sufficiently strong that (3i = 0.06 on average. The 
dominance of the magnetic pressure is again reflected in the narrow, nearly axisymmetric 
sheets of ions surrounded by strong toroidal fleld. At the end of the run, the maximum 
ion density is 8 times the initial value, the largest relative compression of any of this series 
of runs. The ion sheets have a similarity of appearance to the channel solutions seen in 
the single fluid simulations, except that here fleld ampliflcation stops and the growth of 
the nonaxisymmetric parasitic instability is suppressed or signiflcantly reduced. Hence, the 
thin ion sheets endure. Despite a coUisional frequency ~ Q,, the ions and neutrals are not 
sufficient well coupled to overcome the complete dominance of the magnetic fleld on the ion 
evolution. The perturbed ion velocities are larger than those of the neutrals. The neutrals 
are again modifled by the presence of nearly z-independent trailing pressure waves, except 
that the amplitude of these waves is now only SPn/Pn — 0.01. 

The time-averaged Maxwell stress in Z23 is 4.8 x 10~'^{Lz^)'^ , the ion Reynolds stress 
is 5.5 x 10~'^{Lzfl)'^, and the neutral Reynolds stress is 7.5 x 10~^{Lzfiy The stress is 
dominated by the Maxwell component which remains proportional to the magnetic energy. 
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with amag = 0.4. The ion Reynolds stress is proportionally larger than the neutral stress 
by a factor of 7, but, due to the low ion fraction, its absolute value makes no significant 
contribution. 

To summarize, we have examined three decades in ionization fraction with coupling 
frequencies near the orbital frequency. The observed numerical growth rates are consistent 
with the linear analysis of BB. Further, we find that these collision frequencies are 
insufficient to involve the neutrals in the evolution beyond providing a drag. The saturation 
amplitudes are determined by the ion Alfven frequency, and low Pi values. Because 
saturation is determined mainly by the ions, and magnetic equipartition with the ion 
thermal energy, ion-neutral systems with very low ionization fractions and weak coupling 
would be expected to saturate at small magnetic field strengths. 



3.2.3. Effect of the collision frequency 



The analysis of BB highlights the importance of the ratio of the collision frequency to 
the epicyclic frequency. In the simulations described so far we have maintained the collision 
frequency near its critical value, i.e., comparable to the epicyclic frequency. We next 
investigate the effects of altering the collision frequency by modifying the drag coefficient 
7. How does the turbulent magnetic energy vary as the coupling frequency is increased? At 
what point as the coupling frequency is increased do the neutrals couple strongly enough 
with the ions to be significantly affected by the MHD instability? At what value is the 
coupling frequency so low that the ions essentially decouple from the neutrals? 

To investigate these questions, we run simulations with the same ionization fraction 
and initial magnetic field strength, but with different values of the drag coefficient 7. We 
begin with the initial conditions of the fiducial run, / = 0.1 and initial field strength 
VAzi ^ 0.044L^O, and set 7 = 0.001 in run Z24, 7 = 0.1 in Z25, and 7 = 1.0 in Z28. Along 
with the fiducial run, this produces an ensemble of simulations with 70/pi — 0.111, 1.11, 
11.1 and 111.1. Time and space averaged values from these simulations are given in Table 
3. 

The evolution of the magnetic energy in these simulations is shown in Figure 6. Along 
with the two-fluid runs, this flgure includes two single-fluid simulations labeled CI and C2. 
Run CI has p = 0.111 (the same as pi in the flducial run) and C2 has p = 1.111 (the same 
as Pi + pn in the flducial run). We use the flducial run's initial magnetic fleld strength 
for both these simulations, hence CI and C2 have different initial values of the Alfven 
speed, specifically va/^ = 0.044 and 0.139. In both of these control simulations, turbulent 
saturation occurs after three orbits, with time-averaged magnetic energies of 0.025(^2^)^ 
and 0.083{Lzfiy. While these values are somewhat below that predicted from the HGB95 
empirical relationship (6), the ratio of the turbulent magnetic energies in these two runs is 
0.30, consistent with the functional dependence expressed in equation (6). 

The 7 = 0.001 run (Z24) is similar to the low density single fluid calculation (CI). The 
ions are again squeezed into the very narrow filaments between strong magnetic fields while 
the neutrals are barely disturbed. This run is the first case discussed where the ion kinetic 
energy exceeds the neutral kinetic energy. The squeezing of the ions by the magnetic field 
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is sufficiently great that the ions sheets become effectively one zone wide and the simulation 
ends. An isothermal single-fluid simulation behaves like Z24. The CI run uses an adiabatic 
equation of state; the additional pressure keeps the fluid sheets resolved, permitting the 
computation to continue. Again, this demonstrates that compressibility effects become 

important when (3 < 1. 

The differences between the 7 = 0.01 run (Z17) and the 7 = 0.1 run (Z25) are small 
but significant. The magnetic field energies are comparable, but Z25 has a 25% larger 
average radial magnetic field strength, and greater neutral Reynolds stress leading to a 20% 
larger a value. Thus, an increase by 10 in the coupling frequency has had some infiuence 
on the saturation level, but has not yet brought the magnetic energies to near their single 
fluid values. 

When the drag coefficient is increased to 7 = 1 (Z28) the coUisional frequency exceeds 
the orbital frequency by 111. The resulting turbulent magnetic energy is significantly larger 
than in the fiducial run, and continues to grow slowly during the simulation. The average 
values between orbits 10 and 14 are Bl/Sn = 0.0049(L^J1)2, B^/8n = OmO{L^Qy, and 

B^/Stt = 0.0017(L^Q)^, about 70% of the values seen in the single fluid simulation C2. 

Table 2 lists the ratios of the ion to neutral kinetic energies. If the two fluids were 
fully coupled, this ratio, and a similar ratio for the Reynolds stress, would be precisely the 
ionization fraction. This provides one measure of the effective coupling in the turbulent 
state. As the coupling frequency increases, the ratio of the Maxwell stress to neutral 
Reynolds stress increases toward the single-fluid range of 4-5. The ion Reynolds stress also 
increases, but even for run Z28 it is not yet quite one tenth the neutral value. It appears 
that coupling frequencies must be ^ 100 before the evolution becomes similar to that 

of the fully ionized fluid. 

As another test of the role of 7, we compare a series of / = 0.01 runs with 7 values 
equal to 0.01, -fil/pi = 0.1 (Z19), 7 = 0.1, -fil/pi = 1 (Z20), and 7 = 10, 'f^/pi = 100 
(Z27). In the least well-coupled model (Z19) the ions are squeezed into very narrow 
filaments between regions of strong toroidal field. With increasing 7 the ions are less 
confined into coherent sheets, more turbulent, and more like the single-fiuid simulations. 
The angular momentum transport is again dominated by the Maxwell stress. In the three 
runs, ftmag = 0.2, 0.3, and 0.4. The increase in Omag refiects the larger values of radial field 
that are obtained with increasing 7. In Z27, the ratio of the ion to neutral kinetic energy 
and ion to neutral Reynolds stress are nearly equal to the ionization fraction, indicating 
good coupling. In order of increasing 7, the turbulent magnetic energies are 2.6 x 10~^, 
1.0 X 10"^, and 6.7 x 10-^{L^nf. These energies are dominated by the toroidal field; the 
ratios of toroidal to radial magnetic energy are 0.006, 0.009, and 0.027. Although the 
best-coupled run (Z27) has the largest magnetic energy, its value is still a factor of 3 below 
the single-fluid prediction with density ptot- It should be noted that in terms of the total 
gas pressure, the initial fleld is quite weak with /3tot = 10^. For a single fluid simulation with 
this strength field, the critical wavelength is less than the grid spacing, so this simulation is 
somewhat under-resolved. 

Next, we consider a well-coupled model (Z22) with a stronger initial field and a lower 
ionization fraction, specifically, / = 0.001 and 7 = 100. This raises the collision frequency 
to 'jpi/fl = 100. The initial field is set to f3i = 0.1. The instability develops rapidly as 
a channel solution (two channels in the vertical direction), then saturates at 4.5 orbits as 
the channels break up into turbulence. Despite beginning with a strong field relative to 
the ions, the coupling is sufficient to raise the effective inertia of the system, and to ensure 
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that the neutrals become fully turbulent as well. The magnetic energy at the end of the 
simulation is 0.036{Lz^)'^ which corresponds to /5 = 28, and f3i = 0.002. As always, the 
Maxwell stress is the largest with a value of 0.018{Lzfiy, but the neutral Reynolds stress is 
now also significant at 0.0073(1/2^)^. The ion Reynolds stress is smaller by a factor of /. 
The value of a is 0.025, while ctmag = 0.7. Overall the evolution is much like that for a single 
fluid. If the collision frequency is sufficiently large, full MHD turbulence and significant 
angular momentum transport are recovered, even for low ionization fractions. 

To explore the other extreme, and to follow the approach to decoupling in the low-7 
limit, we run 4 models with / = 0.01 and 7 = 0.1, 0.05, 0.01, and 0.001 (Z4a, Z5a, Z6a, 
Z7a). These models use an adiabatic equation of state with an initial ratio of ion to neutral 
gas pressure of 4 x 10~^. The adiabatic equation of state prevents the magnetic field from 
squeezing the ions into numerically unresolvable sheets. The runs of this group are all very 
similar except for Z7a. Although all four runs end up with similar toroidal field energies, 
only Z7a has a substantial radial field energy; this results in a larger total Maxwell stress. 
The average energy in the radial ion velocity fiuctuations is almost two orders of magnitude 
larger in Z7a compared with the others. The observed growth rate in Z7a is comparable to 
the single fiuid case. These results indicate that when 7Pi/f2 < 0.01 the ions and neutrals 
behave essentially as two distinct fiuids. 



3.24. Other effects 



BB determined that in two-fluid systems the presence of a strong toroidal fleld can 
affect growth rates. We demonstrate this point by rerunning the fiducial run with the 
addition of a strong toroidal field with By — lOS^ (YZl). Such a field has VA,p > l-6c^j, and 
with '~fPi/Q = 1.11 the growth rates should be significantly affected (BB). Indeed, we find 
that the growth rate of the perturbed magnetic field drops to 0.18f2. After saturation at 
orbit 7 the radial magnetic and kinetic energies are reduced compared to the fiducial run, 
and this reduces the total transport to a = 0.002. The total magnetic energy is slightly 
larger than in the fiducial run, although this is mainly because here the toroidal field began 
at this level. 

The results discussed so far point to a significant role for compressibihty in the growth 
and saturation of the instability. The ions are strongly affected by magnetic pressure when 
Pi < 1 and the ion- neutral coupling frequency is near or below ~ Q. To further study the 
role of compressibility, we rerun the fiducial case with the ion and neutral masses equal (run 
Z18); this increases the ion pressure. Both this run and the fiducial run arc the same during 
the linear growth phase. Run Z18, however, saturates at a magnetic field energy level that 
is about twice that of Z17. Since the saturated /3i value is larger in Z18 (even though the 
magnetic field is stronger than in Z17), the ions are less tightly confined. This is consistent 
with the conclusion that if the saturated magnetic energies yield /? > 1, the gas pressure 
remains mostly unimportant, whereas if /3 < 1, pressure can bring field amplification to 
an end. The effects of pressure are particularly noticeable in the two-fluid simulations 
because, in contrast to previously studied single-fluid models, the magnetic pressures in 
the saturated states generally have Pi < 1. (This depends, of course, on how we chose to 
initialize the simulations. Here, the initial ion pressure is generally much less than the 
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neutral pressure, and the neutral pressure is chosen so that Cs„ = L^Q.) If the coupling is 
sufficiently weak such that the effective inertia for the magnetic instability is close to that 
provided by the ions alone, the pressure prevents the magnetic field from becoming much 
stronger than equipartition with the ion pressure. 



The structure associated with the nonlinear regime of the instability in the weakly- 
coupled pure Bz field simulations described above is dominated by the persistence of the 
channel solution in the form of narrow, nearly axisymmetric sheets of ions. This is in 
contrast to the single-fluid simulations, where parasitic instabilities inevitably cause the 
channel solution to break up into MHD turbulence (HGB95). The axisymmetric coherence 
of the channel solution is lost when the initial field is not uniform. To study this in 
the two-fluid case, we have performed simulations which begin with fields of the form 
Bz oc sin(a;) so that there is zero net flux through the computational domain. 

Our first zero net B^ field simulation (ZNl) is computed using nearly the same 
parameter values as the fiducial pure B^ field simulation, i.e., /3j = 80, / = 0.1, and 
7 = 0.009, yielding a coupling frequency of ^Pi/Q = 1. We adopt an isothermal equation 
of state, and use the standard resolution to evolve the model. Once again, we observe 
exponential growth of the magnetic energy, with saturation at an amplitude of 0.01(^2^2)^ 
occurring near 5 orbits. In the saturated state, the Maxwell stress dominates the Reynolds 
stress, with —{B^By/An) — 0.002{Lzfiy averaged over the ffist 10 orbits. These values are 
all consistent with those reported for the fiducial run Z17. As is the case in single-fluid zero 
net Bz field simulations (Hawlcy, Gammic, & Balbus 1996), the amplitude of the turbulence 
decreases significantly after saturation, but it never dies completely away; the magnetic 
energy remains at least a factor of 5 larger than in the initial state. 

We have also computed the evolution of a zero net Bz field with a smaller drag 
coefficient, i.e. 7 = 0.002 yielding ^pi/fl = 2/9, but all other parameters identical to ZNl. 
This model (ZN2) evolves in a similar fashion, but with a smaller initial growth rate that 
produces saturation at a slightly later time (6 orbits). The amplitude of the saturated 
magnetic energy, kinetic energy. Maxwell and Reynolds stress are all similar to the pure Bz 
field simulation Z24. However, unlike Z24, the ions do not have a channel-like structure. 

The turbulence observed in these simulations is driven by the instability acting upon 
the ions. Because of drag, turbulent motions are also driven in the neutrals. However, the 
ions and neutrals will be coupled only on scales greater than Ldrag — ^An/lPi- Since in 
these models the magnetic field saturates at vai ~ Csi we can write 



For the parameter values adopted for run ZNl, L^j-^g — 0.09, and for ZN2 L^j-g^g — 0.42. 
To investigate whether there is any quantifiable difference between the turbulence in the 
ions and neutrals on scales above and below Ldrag, we plot in Figure 7 the power spectrum 



3.3. Zero net B^ fields 
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of the specific kinetic energy in the fluctuations of the z-component of velocity {Sv'^) in 
both the ions and neutrals for ZNl (top panel) and ZN2 (bottom panel). In run ZNl, the 
energy in velocity fluctuations in the ions and neutrals is virtually identical on all scales. 
Given the wavenumber associated with the drag length is k^rag = 227r/L2, which is more 
than the largest wavenumber representable at the standard resolution, this result is to be 
expected. However, the power spectrum for run ZN2 shows a systematic difference between 
the ions and neutrals above a wavenumber of about 10; the energy associated with velocity 
fluctuations in the neutrals is less than that of the ions by about an order of magnitude. In 
this case, kdrag = 4:.8n/Lz, which is in agreement with the observed location of the break 
in the power spectrum. In summary, we find that on scales less then L^rag the ions and 
neutrals are poorly coupled, and the structure of the neutrals is very smooth, whereas on 
scales larger than L^rag the ions and neutrals are well coupled and exhibit a decreasing 
power law spectrum with identical amplitude. 



3.4. Uniform By fields 



The analysis of BB considered axisymmctric, vertical wavenumber modes. A remaining 
question is whether purely toroidal fields are also unstable in the two-fiuid system. Toroidal 
fields exhibit a nonaxisymmetric instability in the single-fiuid system, as found both 
analytically and numerically (Balbus & Hawlcy 1992; HGB95; Matsumoto & Tajima 1995). 
In the fully ionized disk, the observed growth rates are lower than with vertical fields, but 
not dramatically so. The turbulent magnetic energy densities that result are lower than 
those seen in the vertical field simulations by about an order of magnitude, although they 
are still proportional to the geometric average of the background field strength and the field 
strength corresponding to the largest possible unstable wavelengths given the dimensions of 
the simulation box, Ly (HGB95). Outward angular momentum transport is produced with 
a total stress that remains proportional to the magnetic pressure, with ctmag ~ 0.5. 

The nonaxisymmetric instability must be present for the two-fiuid system as well, at 
least in completely coupled and decoupled limits. Here we investigate what happens near 
the critical value of the collision frequency, 'jPi/Q = 1 with the scries of toroidal field runs 
listed in Table 4. Figure 8 shows the time evolution of the toroidal and radial magnetic field 
energies for the models listed in Table 4. Table 5 lists the value of some averaged quantities 
after saturation for these runs. 

The first simulation (Y5) has an ion fraction oi f — 0.1 and a drag coefficient 7 = 0.01. 
Again we use an isothermal equation of state, and assume a ratio of ion to neutral mass 
mi/rrin of 30/2.33. The initial field is a uniform By field with strength /3j = 10, which 
corresponds to a value of kvAi/^ — 0.125Ly/X. The perturbations grow slowly in this 
simulation, leveling out after 20 orbits, continuing to grow very slowly after that time. 
Relatively little total field amplification occurs, with the average magnetic field strength 
near the end of the run (46 orbits) equal to 2.1 x 10~^. Total stress is also small, with 
a = 4.4 X 10~^ and ctmag = 0.22. The ions are found in small scale features that have a 
mostly m = 1 structure. The perturbation in the neutrals consists of trailing m = 1 waves 
that are nearly independent of height z. 
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In these simulations, the fastest growing unstable wavelengths are short compared to 
the computational domain size, although the full range of unstable wavelengths is large 
and includes lengths comparable to Ly. Because of this, numerical resolution will play an 
important role in the linear growth and possibly in the final saturation of the instability. 
To measure this effect, we repeat the fiducial toroidal field run with twice the resolution, 
i.e., 63 X 127 x 63 (LYl). The observed growth rate during the initial orbits is increased. 
At the end of the high resolution run the magnetic energy is twice the value at the same 
time in Y5. At this point the growth rate levels off to a value much like the lower resolution 
simulation. Given the low growth rates it seems impractical to carry out this simulation to 
late times, but it is clear that resolution significantly affects the initial evolution of these 
toroidal field models. In analyzing these runs, therefore, we will focus on qualitative effects 
and relative values of quantities. 

In the next model (Y6) we increase the initial By field to /3j = 0.8, which corresponds 
to kvAi/^ — OAAlLy/X. This is the same ion Alfven speed as in the toroidal plus run 
discussed above (YZl). The initial growth rate is larger than in Y5, but at saturation, 
which occurs at orbit 15, the radial field energy has risen only to 1.5 x 10^'^[Lz^iy . Again 
the stress levels are low, with a = 8.4 x lO^'^ and Omag = 0.072. The low stress levels, and, 
in particular, the low level of ftmag is due to the weak amplification of the poloidal field, and 
its small value compared to the initial toroidal field. 

Generally, larger growth rates are obtained when kvj^/Q ~ 1, but with Y6 we are 
already in the regime where the initial field is superthermal (in the ions) . The linear analysis 
of BB and the single fluid analysis of Balbus & Hawley (1992) found that compressibility 
becomes important for such strong fields, reducing the effective growth rate. In the present 
simulations the ion pressure can be increased and the value of (3i decreased for a given 
field strength if the ion and neutral masses are set equal, = m„. We repeat the above 
two runs with this change, carrying out runs with Vj^;^i/Q = 0.441 (Y7; compare with Y6) 
and vazi/^ = 0.125 (Y8; compare with Y5). The increase in ion pressure has relatively 
little effect on the evolution during the linear phase. After 20 orbits, however, the field in 
run Y8 grows to higher levels, suggesting that the ion pressure plays a significant role in 
the nonlinear evolution. This is consistent with the vertical field results discussed above. 
In contrast, comparing runs Y6 and Y7 (with (3i = 0.8 and 10.3 respectively) produces a 
noticeable difference during the linear phase. The growth rate is increased and the total 
field is amplified by a factor of 2.3 to f3i = 4.4. Angular momentum transport is similarly 
enhanced, with a = 0.008 and Oimag = 0.33. These runs confirm that compressibility can 
affect both the linear and nonlinear evolution of the instability whenever /3j < 1 . 

As the collision frequency increases, the system should approach that of a single fluid. 
To test this we repeat the fiducial run with the collision parameter increased to ■jPi/Q =111 
(Y9). In this run the initial growth through 15 orbits is quite similar to the first run (Y5), 
but poloidal field growth continues for a longer time than in Y5. At 40 orbits the field 
energies are still rising, and have twice the magnetic energy as Y5. The tight coupling 
means that the perturbed velocities in the two fluids have the same values on average. 



4. Discussion 
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Wc have carried out a series of simulations to examine the behavior of the 
magnetorotational instabihty in partially-coupled ion-neutral systems. We find that weak 
vertical fields are unstable with growth rates consistent with the linear analysis of BB. The 
initial evolution of the vertical field instability behaves much like the channel solution found 
in the single fluid simulations. This is consistent with the flnding of BB that the unstable 
two-fluid mode with vertical magnetic fleld is still an exact local nonlinear solution to the 
equations of motion. In our simulations, if the ion fraction and the collision frequency 
are low the channel solution can persist, squeezing the ions into narrow channels between 
regions of strong toroidal fleld. As Pi decreases, the simulation ends when the channels 
become one grid zone thick. 

Brandenburg & Zweibel (1995) have emphasized that because ion-neutral coupling 
acts as a nonlinear diffusion process, certain initial field configurations can evolve into thin 
current sheets. Because we are far from the well-coupled limit studied by Brandeburg 
& Zweibel, it is difficult to asses the relative importance of nonlinear diffusion due to 
ambipolar diffusion in these simulations. We have identified the primary mechanism for the 
growth of the narrow sheets in the vertical field simulations with the "channel solution" of 
the magnetorotational instability. Precisely the same thin sheets are observed in single-fluid 
(fully ionized) simulations, and in the ions in the nearly-decoupled limit. In both cases, the 
fluid is squeezed by an exponentially growing fleld until the channel becomes one zone thick 
and the simulation ends. Further, the zero-net-fleld two-fluid simulations presented here, 
which do not produce a channel solution, also avoid the formation of large-scale current 
sheets. It remains an open question whether the nonlinear diffusion of Brandenburg Sz 
Zweibel (1995) increases the tendency to thin sheets, or otherwise affects the turbulent 
state in ion-neutral systems. 

Through a series of simulations with varying drag coefficient 7 we have explored the 
transition from the fully coupled to the completely uncoupled limit. We find that there 
are important two-fluid effects over the range 0.01 < < 100. In simulations at the 

upper end of this range the ions and neutrals are strongly coupled, although the turbulence 
levels remain somewhat below the single-fluid case. Our results are consistent with the 
simulations of Mac Low et al. (1995) who found signiflcant ambipolar diffusion when 
7Pj/n = 30, but little apparent diffusion at a value of 300. 

At the other extreme, the weakly coupled limit with 7Pj/f2 ~ 0.01, the perturbed 
ion kinetic energy actually exceeds that of the neutrals (despite an ionization fraction of 
0.01), and the magnetic energy is increased over levels found in the results for larger drag 
coefficient values. This signals the almost complete separation of the ion and neutral fluid. 

When the coUisional frequency is comparable to the orbital frequency, 7Pi/f^ ~ 1, the 
nonlinear evolution of the instability is primarily controlled by the ion density. Magnetic 
fleld saturation occurs near or just above equipartition with the ion pressure. Hence the 

saturated magnetic fleld energy goes down with ionization fraction, which, in turn, means 
that the total angular momentum transport is also proportional to /. 

The instability transports angular momentum in the two-fluid system, through a 
combination of Maxwell and Reynolds stresses. The total stress is proportional to the 
magnetic energy density in all simulations. The constant of proportionality (ctmag) depends 
on the amount of radial fleld ampliflcation relative to the toroidal fleld. With good coupling 
or larger ionization fractions ctmag is comparable with the single-fluid value (HGB95). Weak 
coupling (but not completely uncoupled) and low ionization fractions produce reduced 
values. 
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Because the neutrals are coupled to the magnetic field only through drag with the ions, 
the efi^ect of the instability on the neutrals depends on the size of the drag coefficient. The 
neutral system alone is hydrodynamic, and as such is stable to hydrodynamic perturbations 
(BHS96). The ion velocity fiuctuations are correlated so as to produce a net positive ion 
Reynolds stress, and that correlation is conveyed to the neutrals. However, because the 
neutrals lack direct coupling to the magnetic fields, a positive Reynolds stress acts as a sink 
in the neutral angular velocity fluctuation equation. In the intermediate coupling regime 
this apparently works to damp the strength of the ion turbulence to levels below what 
would be seen with the ion system fully decoupled. 

Why does gas pressure matter in these simulations and not in the single-fluid studies 
done previously? In the single-fluid simulations (HGB95) the gas pressure was largely 
irrelevant, in part because the box size was chosen so that P = (L^fi)^ and the largest 
unstable wavelength corresponded to a subthermal fleld. In the ion-neutral system as we 
have set it up here, fields that are superthermal in the ion pressure (i.e., have Pi < 1) 
can still be unstable. For small ionization fractions these fields will be very subthermal in 
terms of the total (neutral plus ion) pressure (/3 ^ 1). When j^l/pi ~ 1, the ion pressure 
significantly affects the nonfinear stage of the evolution, for both toroidal and vertical 
initial field models. The tendency for the ions to form thin sheets between regions of 
strong toroidal field is symptomatic. In some cases this reduces the vigor of the resulting 
turbulence. Stronger turbulence and more angular momentum transport can be obtained 
by decreasing the ion mass, and hence increasing the ion pressure (under the assumption of 
equipartition of thermal energy between the ion and neutral particles, mjC^j = m„c^j). We 
consistently find stronger turbulence with higher ion pressures. 

When toroidal fields are included in the initial conditions along with vertical field, 
they can reduce linear growth rates, in accordance with the analysis of BB. However, the 
reduced growth rates have less of an effect once the nonlinear regime is reached. In fact, far 
from preventing the the growth of the instability, we find that purely toroidal initial fields 
are also unstable. However, the instability is significantly reduced for weak coupling if the 
toroidal field is superthermal. Although this is true in the single fiuid case as well, there 
it was hardly an issue because va ~ Cg was already a very strong field. Here, however, if 
rrii > rrin, the ion pressure can be substantially smaller than the neutral. If the ion fraction 
is small, the ion Alfven speed can easily become large compared to the ion sound speed, 
despite a relatively weak toroidal field. 

Numerical resolution is more important in these simulations than for the single fluid 
models. The presence of both ion and neutral components introduces additional Icngthscales 
that can be quite disparate, as in the case of low ionization fractions and intermediate 
strength coupling. A Fourier analysis of the velocity fluctuations associated with the ions 
and neutrals shows that, as expected, on scales greater than L^rag = VAn/jPi, the ions and 
neutrals are well coupled with identical power spectra for velocity fluctuation, whereas on 
scales less than L^h-ag the fluctuations in the neutrals are at a considerably lower amplitude 
than those in the ions. Thus, Ldrag represents an important lengthscale introduced into two 
fluid models that must be resolved numerically. 

In summary, we expect that the two-fluid effects studied here will be appropriate to 
a transition region in protostellar or protoplanetary accretion disks between the inner, 
hot, and fully ionized regions, and the outer, cold, and essentially neutral regions. Our 
simulations suggest a more stringent criterion for good coupling than that obtained from 
the linear result of BB. We flnd that, while the linear instability is present for coupling 
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frequencies ~ fl, significant turbulence and angular momentum transport can only occur 
when coupling frequencies jfl/pi > 100. For weaker coupling, the magnetic energies are 
determined primarily by the ion density. 

There are many interesting issues regarding this transition region in a protostcllar 
disk that go beyond the limitations of this initial study. Here we have assumed that both 
the recombination and ionization timescales are much longer than the orbital period. In 
real systems, the ion density in structures such as the ion filaments may be limited by 
recombination. On the other hand, at the interface between weak and strong coupling, 
increasing the ionization fraction increases the efficiency of magnetic coupling, which in 
turn increases the level of turbulent heating thus raising the ionizational level yet higher. 
Investigating stability issues such as these, as well as the global structure and dynamics of 
weakly coupled disks, are fruitful areas for future research. 

We thank Steven Balbus for valuable discussions, and Omer Blaes and Mordccai-Mark 
Mac Low for comments on the manuscript. This work is supported in part by NASA 
grants NAG-53058, NAGW-4431, and NSF grant AST-9423187 to J.H., by NASA grant 
NAG-54278 and NSF grant AST-9528299 to J.S., and by a metacenter grant MCA95C003 
from the Pittsburgh Supercomputing Center and NCSA. 
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TABLE 1: Z FIELD SIMULATIONS 



Model 


f 


0' 




m,- /m„ 






Grid 




Orbits 


Comment 


Z17 


0.1 


0.01 


1.1 


12.9 


0.044 


31 X 63 X 31 


26 


Fiducial 


LZl 


0.1 


0.01 


1.1 


12.9 


0.044 


63 X 127 X 63 


13 


High- Res 


Z20 


0.01 


0.1 


1.0 


12.9 


0.044 


31 


X 


63 


X 


31 


30 




Z21 


0.001 


1 


1.0 


12.9 


0.044 


31 


X 


63 


X 


31 


71 




Z23 


0.001 


1 


1.0 


12.9 


0.46 


31 




63 




31 


10 




Z24 


0.1 


0.001 


0.11 


12.9 


0.044 


31 


X 


63 


X 


31 


4 




Z25 


0.1 


0.1 


11.1 


12.9 


0.044 


31 


X 


63 


X 


31 


11 




Z28 


0.1 


1 


111 


12.9 


0.044 


31 


X 


63 


X 


31 


14 




Z19 


0.01 


0.1 


0.1 


12.9 


0.044 


31 


X 


63 


X 


31 


27 




Z27 


0.01 


10. 


101 


12.9 


0.044 


31 


X 


63 


X 


31 


29 




Z22 


0.001 


100 


100 


12.9 


1.25 


31 


X 


63 


X 


31 


63 




Z4a 


0.01 


0.1 


1.0 


2.1 


0.044 


31 


X 


63 


X 


31 


32 


Adiabatic 


Z5a 


0.01 


0.05 


0.5 


2.1 


0.044 


31 


X 


63 


X 


31 


32 


Adiabatic 


Z6a 


0.01 


0.01 


0.1 


2.1 


0.044 


31 


X 


63 


X 


31 


24 


Adiabatic 


Z7a 


0.01 


0.001 


0.01 


2.1 


0.044 


31 


X 


63 


X 


31 


19 


Adiabatic 


YZl 


0.1 


1 


1.1 


12.9 


0.044 


31 


X 


63 


X 


31 


14 




Z18 


0.1 


0.01 


1.1 


1 


0.044 


31 


X 


63 


X 


31 


15 
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TABLE 2: TIME- AND VOLUME- AVERAGE VALUES 
EFFECT OF IONIZATION FRACTION 



((Quantity)) 


Z17 




Z20 




Z21 




Z23 




/ 


0.1 




0.01 




0.001 




0.001 




S^/87r(L,0)2 


5.3 X 10" 


-4 


7.8 X 10- 


-6 


1.8 X 10- 


-7 


5.7 X 10- 


-5 




1.1 X 10- 


-2 


1.1 X 10- 


-3 


2.0 X 10- 


-4 


1.1 X 10- 


-3 




3.8 X 10" 


-4 


1.4 X 10- 


-5 


1.2 X 10" 


-6 


1.2 X 10- 


-4 


{pviy/2{L,nf 


6.8 X 10" 


-4 


4.7 X 10- 


-6 


7.0 X 10" 


-7 


8.1 X 10- 


-7 




3.0 X 10- 


-4 


5.5 X 10- 


-7 


7.7 X 10- 


-8 


1.4 X 10- 


-7 




7.0 X 10- 


-5 


3.2 X 10- 


-7 


8.1 X 10- 


-8 


1.3 X 10- 


-7 




0.10 




0.010 




0.0010 




0.0031 




{p5vl)i/ {p5vl)^ 


0.16 




0.013 




0.0010 




0.0048 




{pvl)-J{pvl)^ 


0.13 




0.013 




0.0017 




0.0026 




Em/KE, 


11.8 




200 




5100 




1200 




-B,By/A7T{L,nf 


3.4 X 10- 


-3 


1.4 X 10- 


-4 


9.0 X 10" 


-6 


4.8 X 10- 


-4 


{pvxVy)i/{L^rtf 


3.5 X 10- 


-4 


1.4 X 10- 


-6 


1.8 X 10" 


-8 


5.5 X 10- 


-7 


{pVxVy)^/{L^Vtf 


2.4 X 10- 


-3 


1.2 X 10- 


-4 


1.7 X 10- 


-5 


7.5 X 10- 


-5 


Max/ Reynoldsi ^ 


9.7 




100 




510 




870 




Max/ReynoldSn 


1.4 




1.3 




0.53 




6.4 




a 


6.1 X 10- 


-3 


2.6 X 10- 


-4 


2.6 X 10- 


-5 


5.6 X 10- 


-4 


Q^mag 


0.61 




0.23 




0.13 




0.42 





"Max/Reynoldsi is the ratio of the Maxwell stress to ion Reynolds stress. Similarly 
Max/ReynoldSn is the ratio of the Maxwell stress to neutral Reynolds stress. 
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TABLE 3: TIME- AND VOLUME- AVERAGE VALUES 
EFFECT OF COUPLING CONSTANT 



((Quantity)) 


Z24 




Z17 




Z25 




Z28 




7 


0.001 




0.01 




0.1 




1.0 




Bl/STiiL.nY 


2.4 X 10" 


-3 


5.3 X 10" 


-4 


8.1 X 10" 


-4 


4.9 X 10- 


-3 


y' \-^z""/ 


1.9 X 10" 


-2 


1.1 X 10- 


-2 


1.1 X 10- 


-2 


5.0 X 10- 


-2 


Z 1 \ Z / 


3.7 X 10" 


-4 


3.8 X 10- 


-4 


3.7 X 10- 


-4 


1.7 X 10- 


-3 


\r xJ^l V -2 y 


1.2 X IQ- 


-3 


6.8 X IQ- 


-4 


5.1 X 10" 


-4 


1.3 X IQ- 


-3 




1.7 X 10" 


-3 


3.0 X 10" 


-4 


1.4 X 10- 


-4 


5.8 X 10- 


-4 


{fyvl)J2iL,Vtf 


2.3 X 10" 


-4 


7.0 X 10- 


-5 


8.7 X 10- 


-5 


4.1 X 10- 


-4 




1.04 




0.10 




0.11 




0.11 




[P'J%)i/[P'J%)n 


0.41 




U. iu 




n 1 

U. lO 




n 1 9 




{pv%/{pv% 


0.65 




0.13 




0.12 




0.11 




Em/KE, 


6.1 




11.8 




18 




26 






1.2 X IQ- 


-2 


3.4 X 10" 


-3 


4.7 X 10" 


-3 


2.4 X 10- 


-2 




2.2 X 10" 


-3 


3.5 X 10" 


-4 


2.3 X 10- 


-4 


6.4 X 10- 


-4 


{pVxVy)n/{L^QY 


7.6 X 10- 


-5 


2.4 X 10- 


-3 


1.8 X 10- 


-3 


5.5 X 10- 


-3 


Max/Reynoldsi 


5.2 




9.7 




20 




37 




Max/ Reynoldsn 


150 




1.4 




2.6 




4.3 




a 


1.4 X 10" 


-2 


6.1 X 10- 


-3 


6.6 X 10- 


-3 


3.0 X 10- 


-2 


Q^mag 


0.63 




0.61 




0.54 




0.53 
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TABLE 4: Y FIELD SIMULATIONS 



Model 


/ 


7 




rrii/mn 


VAzi/ 


Grid 


Orbits 


Y5 


0.1 


0.01 


1.1 


12.9 


0.125 


31 X 63 X 31 


46 


LYl 


0.1 


0.01 


1.1 


12.9 


0.125 


63 X 127 X 63 


19 


Y6 


0.1 


0.01 


1.1 


12.9 


0.441 


31 X 63 X 31 


32 


Y7 


0.1 


0.01 


1.1 


1 


0.441 


31 X 63 X 31 


17 


Y8 


0.1 


0.01 


1.1 


1 


0.125 


31 X 63 X 31 


50 


Y9 


0.1 


1 


111 


12.9 


0.125 


31 X 63 X 31 


42 
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TABLE 5: TIME- AND VOLUME- AVERAGE VALUES 
Y FIELD SIMULATIONS 





Y5 




Y6 




Y7 




Y8 




'\7'n 
Y9 






1.3 X 10' 


_5 


1.5 X 10" 


-4 


1.1 X 10' 


_3 


7.1 X 10' 


-5 


2.1 X 10' 


-4 




2.1 X 10" 


_3 


1.2 X 10" 


_2 


2.4 X 10" 


_2 


4.2 X 10' 


_3 


5.6 X 10" 


_3 




2.6 X 10" 


-g 


8.9 X 10" 


-5 


4.0 X 10" 


-4 


1.7 X 10" 


-5 


6.8 X 10" 


-5 


( 2\ lc\{ T r\\1 


5.9 X 10" 


-5 


5.9 X 10" 


-5 


5.6 X 10" 


-4 


1.7 X 10" 


-4 


1.5 X 10" 


-4 




1.6 X 10 


-6 


3.8 X 10 


-5 


3.3 X 10 


-4 


3.7 X 10 


-5 


4.1 X 10 


-5 




1.7 X 10 


-6 


1.7 X 10 


-5 


1.1 X 10 


-4 


1.2 X 10 


-5 


3.1 X 10 


-5 


/ 9 \ / / 9 \ 
{P%)i/{pvi)n 


0.094 




0.097 




0.10 




0.11 




0.11 




/ c"9\ // r9\ 


0.11 




0.19 




0.21 




0.15 




0.11 




{pvl)-J{pv% 


0.12 




0.10 




0.096 




0.11 




0.11 




Em/KE, 


35 




113. 




27. 




25 




27. 




-B,By/A7r{L,nY 


1.8 X 10" 


-4 


5.6 X 10" 


-4 


5.7 X 10" 


-3 


7.6 X 10" 


-4 


1.4 X 10" 


-3 


{pV:,Vy)i/{L^nY 


2.4 X 10" 


-5 


3.4 X 10" 


-5 


3.3 X 10" 


-4 


5.8 X 10' 


-5 


5.5 X 10' 


-4 


{pv^Vy)J{L^Q.Y 


2.4 X 10" 


-4 


2.5 X 10" 


-4 


2.4 X 10' 


-3 


4.5 X 10' 


-4 


5.0 X 10' 


-4 


Max/Reynoldsi 


7.6 




17. 




17. 




13. 




26. 




Max/ReynoldSn 


0.75 




2.3 




2.3 




1.7 




2.9 




a 


4.4 X 10" 


-4 


8.4 X 10" 


-4 


7.6 X 10' 


-3 


1.1 X 10' 


-3 


2.0 X 10' 


-3 


'-''mag 


0.22 




0.07 




0.33 




0.30 




0.34 
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Fig. 1. — Numerical growth rates for the radial field in two-dimensional simulations with 
ionization fraction / = 0.1. The sohd squares are the 64 x 64 grid zone simulations (5 cases); 
the open stars are the 32 x 32 grid zone simulations (three cases). The sohd line is the 
analytic growth rate from linear stability theory. 

Fig. 2. — Time evolution of the individual components of the magnetic energy for the 
fiducial run Z17 (bold lines) and the high resolution version of the fiducial run (LZl). 

Fig. 3. — Volumetric rendering of (a) the ion density and (b) the magnitude of the toroidal 
field in the fiducial run at orbit 15. In (a) brightness is a function of density, whereas in (b) 
the dark regions correspond to strong field. Comparing the two figures shows that the ions 
lie in thin sheets sandwiched between regions of strong field. 

Fig. 4. — Time evolution of the ion and neutral perturbed kinetic energies in the fiducial 
run (Z17). The behavior of the two curves is very similar; they are offset by the neutral/ion 
density ratio. 

Fig. 5. — Time evolution of the magnetic energy, normalized as the Alfven speed squared for 
runs Z17, Z20 and Z21 which have ionization fractions of 0.1, 0.01, and 0.001, respectively. 
While the total toroidal field amplification is comparable in the three runs, lower ionization 
fractions produce smaller poloidal field amplification. 

Fig. 6. — Time evolution of total magnetic energy for a series of runs with increasing drag 
coefficient 7. The curves are labeled by their run number, as hsted in Table 1. Also included 
are two single-fiuid comparison runs, CI (with density corresponding to the ion density in 
the fiducial run) and C2 (with density corresponding to the total density in the fiducial run). 
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Fig. 7. — Power in fluctuations in Vzi (solid) and Vzn (dashed) versus wavenumber in the 
y-direction for the zero net field runs ZNl (top panel) and ZN2 (bottom panel). The "drag 
length" is given by eq. (7); L^rag — 0.09Lz for the top panel, and L^rag — 0A2Lz for the 
bottom. In ZN2 there is more power at high wavenumbers in the ions than in the neutrals, 
while in ZNl they are both comparable. This agrees with the expectation that when the 
drag length is large (bottom), the ions and neutrals are only weakly coupled, and the ions 
should show more small scale structure. 



Fig. 8. — Time evolution of toroidal (top) and radial (bottom) magnetic field energies in 
the initial toroidal field runs. The curves are labeled by the run number as listed in Table 4. 
All show field amplification, although at lower growth rates than the vertical field models. 
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